/* Erathosthenes.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */


#ifndef Eratosthenes
#define Eratosthenes
#include <iostream>
#include <set>
#include <list>
#include <src/miscellaneous/io/IOFacilities.h>

/*! Eratostene sieve algorithm for a totally ordered ring ( +,-,*,/,%,<).
	This is a generic implementation of the sieve algorithm : 
	consider a list of the first N number, 1 being the first prime number.
	Consider the second prime number 2. Erase from the list all numbers i such that i%2 = 0.
	The next in the list (designing it with 3) is a prime number. Erase from the remaining list all the number such that i%3 = 0, and so on..
	
    There is an hidden difficulty involving performance issues : where should we start the sieve algorithm ?
	For integer values, given a prime number i, we can start at the first remaining number greater or equal to i*i. But this is no longer true
	for general ring, because the order relation i<j does not imply i*k < j*k. For instance on mod 2 polynomial ring
	we have (x+1)*(x+1) = x^2+1 < x*(x+1) = x^2+x, but x+1 > x. This is due to the fact that the addition does not automatically preserve ordering (x+x = 0).

	However if i%j = 0, then i*k < j*k. Indeed in this case we have trivially (i*k) % (j * k) = 0, i.e. i*k < j*k (is this really true ?).
	Thus, if i is the next prime, we start the sieve at the first number j verifying j%i = 0.
	
*/
template <class K>
	class Erathostenes 
	{
		public :
			typedef K										Body;
			typedef typename std::set<K>::iterator			iterator;
			typedef typename std::set<K>::const_iterator	const_iterator;
//			Erathostenes(): first_(1), zero_(0){};
			Erathostenes( const Erathostenes & from) : first_(1), zero_(0), first_list_(from.first_list_){
				first_it_ = first_list_.find(*from.first_it_);
			};
			Erathostenes& operator = ( const Erathostenes & from) {
				first_list_ = from.first_list_;
				first_it_ = first_list_.find(*from.first_it_);
				return *this;
			};
//! The first prime are computed. Default is a buffer of N=100 numbers.
			Erathostenes(unsigned long N = 100) : first_(1), zero_(0), N_(N){
				K value(first_);
				for (unsigned long i = 1; i < N_; ++i)
							first_list_.insert(first_list_.end(),value++); // postfix here
				iterator it_first = first_list_.begin();
				++it_first;
				iterator it_from(it_first);
				bool is_not_prime = true;
				for (++it_from; it_from != first_list_.end();) 
				{
					K breakpoint1 = (*it_from);
					K breakpoint2 = ( (*it_first) * (*it_first) );
					K breakpoint3 = breakpoint1 / breakpoint2;
					is_not_prime = breakpoint3 != zero_;

					if (is_not_prime == true)
					{
						erase(it_first,it_from--); //it_from will be deleted in erase.
						++it_first;
					};
					++it_from;
				}
				first_it_ = first_list_.begin();
		};
//! get the next prime number. If we overpass the prime number buffer N, we recompute the prime numbers for N:=2*N.
			const K & operator++() {
				++first_it_;
				if (first_it_ == first_list_.end()) {
					K temp = *first_list_.rbegin();
					*this = Erathostenes<K>(2*N_);
					first_it_ = first_list_.find(temp);
					++first_it_;
				}
				return *first_it_;
			}

			void erase(iterator itl,iterator it_from)
			{
				K value(*itl);
//				std::cout << std::endl << "remainder " << first_list_ << std::endl;

				while ( it_from != first_list_.end() ) {
					if (*it_from % value == zero_)
						first_list_.erase( it_from++ );
					else
						++ it_from;
				};
//				std::cout << std::endl << "remainder " << first_list_ << std::endl;
			};
			typename std::set<K> first_list_;
			const_iterator first_it_;
			const K first_;
			const K zero_;
			const std::size_t N_;
		};
 template <typename K>
	std::ostream& operator<< (std::ostream& o, const Erathostenes<K> & t)
	{
			return o << t.first_list_;
	};


#endif // Eratosthenes

